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Abstract. The crystalline ground state of macroions confined between two neutral 
parallel plates in the presence of their homogeneously spread counterions is calculated 
by lattice-sum minimization of candidate phases involving up to six layers. For 
increasing macroion density, a cascade of solid-solid transitions is found involving 
various multilayered crystals. The cascade includes triangular monolayer and buckled 
bilayer as well as rhombic, squared and triangular phase structures. 
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1. Introduction 

Strong correlations in Coulomb systems lead to a variety of new effects which are absent 
for neutral particles, see e.g. Refs. [U [2] for a review. Among those are nonlinear 
screening effects [31 [H [5l [6] , charge inversion [7j, Coulomb criticality [8l[9j, like-charge 
attraction for multivalent ions [iQl [HI [121 [E] as well as exotic binary crystalline 
structures unknown for uncharged systems [HI [15] . 

By using charged colloidal suspensions [I6] or dust particles in plasmas [IT], it is 
possible to realize strongly asymmetric mixtures of oppositely charged particles. These 
systems consist of mesoscopic highly charged "macroions" and microscopic counterions 
with a low valency resulting in strong charge and size asymmetries. Since the charges of 
the macroions are high, strong Coulomb correlations are typical for macroions. Most of 
the physics can still be encaptured by viewing these systems as strongly asymmetric and 
strongly coupled electrolytes. In recent years, it was possible to confine macroions in 
sheets between two parallel plates [HI [IHl [lEl [20l [21] and to observe the resulting lateral 
structure of the particles. The gross features can be understood in terms of an (effective) 
one-component system with a Yukawa pair interaction [22l [23l [24l [25l [26l [27]. In fact, 
the mono- and bilayer ground-state structures which were obtained from a Yukawa 
model [28| describe the experimentally found structures [29]. For multilayers beyond 
the bilayer regime, a rich variety of stable phases are found in experiments [30l [HU [32] 
as well as in simulations [33] , which are all theoretically confirmed for a Yukawa system 
between two neutral walls [34]. This motivates a study about the influence of the wall- 
particle interaction on the phase behaviour of multilayered crystalline sheets in slit-like 
confinement [35] . 

In this paper, we consider a model for macroions confined between two parallel 
neutral walls There is a direct Coulomb interaction between the point-like particles. 
The total system is charge-neutral and the counterions are kept at high temperature and 
are homogeneously spread between the plates resulting in an attraction acting on the 
macroions towards the middle of the plates. The system is realized for highly charged 
colloidal particles or dust particles in plasmas. Some early theoretical and simulational 
investigations on clusters of artificial atoms [371 [33 [33] and dusty plasmas [iQl [17] 
as well as one-component plasmas [41], [421 [13], including all the parabolic potentials 
acting as confinement, reveals the existence of multilayers. We therefore include the 
regime beyond bilayers in our discussion. Lattice sum minimizations among a broad 
set of candidate structures are used to determine the structure which minimizes the 
potential energy per particle. For increasing macroion density, we find a cascade of 
solid-solid transitions which includes triangular monolayer, buckled bilayer and squared, 
rhombic and triangular bi-, tri-, tetra-, penta- and hexalayers §. Comparing the results 
to those involving a Yukawa interaction [M], we show that the topology of the phase 
diagram depends crucially on the particle- wall interaction. In fact, some complicated 

:|: Different from [36j we include here a neutralizing background of counterions 
§ For colloid-polymer films, see [44] 
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tetralayered structures which were found stable for the confined Yukawa model are 
unstable in the present model. The strong correlation between phase behaviour and 
wall-particle interactions suggests to tailor new crystalline structures (e.g. with desired 
filtering properties [45]) by a suitable surface treatment of the plates. 

The paper is organized as follows: the model is introduced in section II. After 
discussing the structure of different crystalline multilayers, results for the cascade of 
solid-solid transition are presented in section III. Finally we conclude in section IV. 



2. The Model 

We consider classical point-like particles of charge q {macroions) interacting via the 
unscreened Coulomb pair potential 

V{r) = ^ (1) 

er 

where r denotes the interparticle distance and e the (relative) dielectric constant of 
surrounding medium. The system is confined between two parallel hard walls of area 
A and separation L, see figured! The global charge neutrality of the system is ensured 
by counterions. The latter are taken into account by an homogeneous neutralizing 
background that is smeared out over the whole slit. We mention that we neglect the 
discrete nature of the counterions in this approach, as well as any local ion-counterion 
coupling. 
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Figure 1. A schematic illustration of the model. The ions (e.g., charged colloids) are 
represented by filled circles. The counterions are smeared out between the two hard 
walls located at 2; = ±L/2. This charge distribution generates a quadratic potential 
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in between as shown. The separation between outermost layers (dashed 



lines) is denoted by D. 
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As a consequence of Gauss law, the electric field £"5 (stemming from the neutralizing 
background) is linear in z inside the slit and constant outside the slit. More specifically, 
we have 

for -L/2 <z< +L/2, 

^^(^) = <; ^^Nn^. (2) 

else. 

We thereby implicitly neglect image charge effects [16], meaning that we assume that 
there is no dielectric contrast at the interfaces (at z = ±L/2). The resulting electrostatic 
potential $5, verifying the matching condition at z = ±L/2, then reads 

f for -L/2 < z < +L/2, 

llL^"^"-2iL 

where the reduced density 

V = (4) 

was introduced. Hence, the potential of interaction ^(r) between a macroion and the 
counterion background is merely given by 

Vhiz) = q^.iz). (5) 

We are now in a position to write the total potential energy per particle as || 

N N N 
i=l j=l 1=1 

In its appropriate rescaled form, u reads (within the slit) 

N N ^ ^ N 



i=\ j=l U i=l 

with r^j = rij/L and = Zi/L^ showing that at prescribed confinement width L the 
energy of the system depends only on 77. Consequently the phase diagram at zero 
temperature is given as a function of 77. 

At each given density ry, we have performed lattice sum minimizations for a broad 
set of candidates of crystalline lattices. In order to handle the long ranged Coulomb 
potential, we have used the Lekner summation method [47j for three-dimensional 
systems with two-dimensional periodicity [48], see also [49j. More explicitly, we 
consider in this work three-dimensional crystals with two-dimensional periodicity in 
X- and ^-direction whose primitive cell is a parallelepiped containing n particles. This 

II To remedy the divergence occurring with the first term of ([6|), a two-dimensional neutralizing 
background is introduced in the Lekner (or equivalently Ewald) sum. This neutralizing background 
(implicitly present in the Lekner and/or Ewald sum) has to be distinguished from the one that we use 
to model the counterions, which is smeared out over the whole volume of the slit. 
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parallelepiped is spanned by the three lattice vectors a = a(l, 0, 0), b = a7(cos 9^ sin 9^ 0) 
and c = i^(0,0, 1), where 7 is the aspect ratio (7 = |b|/|a| = b/a) and 9 is the angle 
between a and b. Furthermore, the n particles are distributed, not necessarily evenly, 
on m layers in the z-direction such that c = |c| corresponds to the distance between 
outermost layers (see also figure ^ . Hereby we restrict ourselves to layered situations 
with an up-down inversion symmetry in the averaged occupancy refiecting the up- 
down symmetry of the confining slit. Under this sole restriction, we consider possible 
candidates with n = 1, • • • , 8 and m = 1, • • • , 6 up to symmetric six-layer structures 
with a basis of up to 8 particles. Furthermore, we also examine the stability of several 
asymmetric buckling phases, as predicted in [50]. For given ry, the total potential energy 
per particle is minimized with respect to the particle coordinates of the basis and the 
cell geometry (7 and 9). The resulting stability phase diagrams are shown and discussed 
in the following sections. 



3. Mono- and bilayer phase behavior 

3.1. Phase diagram 

An increase of r] within the mono- and bilayer regime reveals the existence of five stable 
crystalline mono- and bilayers: lA (triangular), 3 A (staggered triangular), 20 (square), 
2R (rhombic) and 2 A (staggered triangular). The integers indicate the number of layers. 
For increasing 77, the stability cascade therefore reads: 

lA ^ 3A ^ 20 ^ 2i? ^ 2A. (8) 

Most of these phases, corresponding to Wigner crystals predicted in earlier 
theoretical investigations [36l EE], are also found in experiments on charged colloidal 
suspensions jSU [52] as well as in Monte Carlo simulations of confined hard spheres [53j . 
The detailed phase diagram is reported in figure [2l 
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Figure 2. Stability phase diagram of crystalline mono- and bilayers. The five stable 
phases lA, 3A, 20, 2R and 2A correspond to Wigner crystals, found in earher 
investigations (see text for details). Note that the monolayer-trilayer transition occurs 
at 77 ~ 1.37. 



We emphasize that the 3A phase (staggered in an ABC manner, see also table 
d]) intervenes between lA and 20 rather than a buckled phase which is present in a 
situation where the external potential has a vanishing curvature at the origin. 

At small reduced densities 77, particles tend to stay in the potential minimum 
(cf. figure [H) created by the counterion background. This is precisely the origin of 




Figure 3. Order parameter h in the transition regime 1 A to 20 via 3A. The monolayer 
lA buckles at a critical density rjc ^ 1.360901 to a trilayer. 



the stability of monolayered Wigner crystals, which never occurs in purely unscreened 
Coulomb systems. % The triangular monolayer lA is stable up to rj = 1.37. At 
larger densities the mutual repulsive interparticle interactions, first term in equation 
([71) , dominates the competition between the interparticle (macroion-macroion) repulsion 
and particle-background (macroion-counterion) attraction. 

The structure with triangular base shape 3A appears as the first stable multilayer 
(see figure [3]), interpolating between lA and 20. The associated order parameter, 
namely the reduced separation 

between the mid-plane and the outer macroion layer (see also figure [1]) , is continuous at 
the transition lA 3A but discontinuous across the 3A 20 transition, see figure [3] 
and [53] . 

By further increase of ry, one recovers the rhombic phase 2i?, which is continuously 
achievable from the square phase 20 by changing 6^, as indicated in the inset of figure 
m The two geometrical order parameters h and sin0, see figure [H indicate thereby 
a continuous transition for 2B 20. On the other hand, at larger values of ry, the 
transition 2R 2 A is of first order as signaled by the jumps of the two geometrical order 
parameters h and sin0, see figure [H The staggered triangular phase 2A corresponds to 
the ultimate stable structure in the high density regime of bilayers. 

^ Indeed, we found that a rectangular bilayer with size ratio 7 = \/3, proposed as a stable structure 
for very small rj in [36], is always energetically beaten by a buckled {2B) bilayered phase. Seen from 
the top, this structure corresponds to the triangular lattice. 




Figure 4. Order parameter h in the transition regime 20 to 2A via 2R. The 
discontinuity Ah in the developing of the layer-layer separation by the transition 
2R 2A is also shown for clarity. In the inset one can regard how changes in 
the same regime. Corresponding structures are also sketched in the inset. Different 
colors indicate different layers. 



3.2. From monolayer to trilayer - An analytic approach 

We now would like to address the transition lA 3A analytically. To do so, we apply 
a Taylor expansion to u{h) around h = D/2L = 0^ see the Appendix for details. The 
resulting asymptotic expression for small interlayer distances h reads 

""^^^ Bo^ + B,r,'/'h' + B,v'^'h' + Uvh'. (10) 



gVeL ' " " ' 3 

with 

So = -1.960516..., 5i = -3.590668 ... , ^2 = 4.968827 ... . (11) 

The profile of the reduced half layer-layer distance h{r]) is obtained upon minimizing u 
with respect to h, i.e. du/dh = 0, leading to 

It is now a simple matter to obtain the reduced density rjc at which the monolayer- 
trilayer transition (lA 3A) takes place. The mathematical condition is thereby 
h{r] = r]c) = yielding 

An 

VVc = -77^ ^Vc = 1.360901 . . . , (13) 

OlDl 

which is in quantitative agreement with the lattice sum minimization results from 
previous section, see figure O 
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Figure 5. Plot of equation equation ([T4j) (dashed line) and numerical calculations 
for finite h{r]) (solid line) based on full lattice sum minimization near the monolayer- 
trilayer lA ^ 3 A transition. 



By inserting the expression f fT3l ) of r]c in f fT2l ) one obtains 



hyri) = -A (14) 

Noticing that the last denominator in equation (fT^ can be approximated (valid in the 
relevant limit r] ^ rj^) by 277^, we obtain a square-root singularity: 



hm h{v) = J--r^,{v - Vc)'/' -{V- VcY/'- (15) 
This theoretical prediction (fT^ is visualized in figure [51 

4. Multilayers 

The presence of the neutralizing background allows the formation of multilayers with 
m > 3 for large enough densities 77, which is forbidden in the absence of a background 
+ . The physical origin of the stability of multilayers in the present system at large 77 
is basically a balance between the mutual unscreened macroion-macroion repulsion and 
the attractive macroion-background interaction. 

+ There is a simple and clear electrostatic argument to explain the exclusive stability of bilayers for 
charges confined between (charged or uncharged) hard walls without neutralizing volume background. 
One has to note that two equally charged walls do not generate any electric field within the slit, 
and consequently do not alter the stable structure obtained at any other surface charge (including 
neutral walls). Hence, if one considers the special case of two walls corresponding to two-dimensional 
neutralizing backgrounds where the ground-state is the 2 A bilayer, we deduce from this that the ground 
state structure is always a bilayer. 
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Figure 6. Stability phase diagram of crystalline multilayers in the presence of a 
neutralizing background. SD, 3i?, 3 A, 40, 4i?, 4 A, 5i?, 5 A and 6i? are obtained as 
stable in the analyzed r^-regime. The corresponding structures are given in table [TJ 



We shall now analyze in detail the high density regime up to 130. Beyond the 
bilayer regime, that is limited by 2A, the cascade found here upon increasing rj reads: 

• • • 30 ^ 3i? ^ 3A ^ 4n ^ 4i? ^ 4A ^ 5i? ^ 5A ^ 6i? • • • , (16) 

where rhombic phases 3i?, 4i?, 5R and 6R have the stacking sequence ABA^ ABAB^ 
ABABA and ABABAB while the triangular phases 3 A, 4 A and 5 A occur as ABC^ 
ABC A and ABCAB^ respectively. More structural details are given in table [H The 
corresponding phase diagram is depicted in figure [61 



Table 1: Structural details and schematic illustration of the stable 
crystalline multilayers. The layers are labeled as follows. The 
bottom one located at z = —D/2 corresponds to first layer (labeled 
as i = 1), and the labels of the successive layers are incremented 
accordingly. For m > 3, the separation between the two first 
layers is characterized by 5D with l/(m — 1) < 5 < 0.5. The 
relative separation vector between two particles of a primitive cell 
belonging to two layers i and j is given by dij. For six layers, the 
separation between the first and the third layers is specified by XD 
with 2/5 < A < 0.5. In the top views of 3 A, 4A, 5 A and 3R, 4R, 
5i?, 6R each basis shape (triangular or rhombic) is emphasized with 
white lines. The rhombic stripes of 3i?, 4i?, 5R and 6R are shown 
again in corresponding perspective views, for clarity. Particles from 
different layers are identified by different colors. 
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The primitive cells of all stable phases found in this work consist of one particle 
per layer. Each constitutive layer possesses the same basis shape (A, □ or i?). These 
layers are shifted to each other, see table [H Note that (for m > 3) the layers become 
equidistant only in the limit ry ^ oc. A remarkable finding is the absence of prism 
phases (at m = 4) that are encountered in hard sphere systems [30l |33] and Yukawa 
systems at finite screening [34] . 

A further overview of the full phase diagram ranging from triangular monolayer 
to rhombic hexalayer structures is shown in figure [71 where the profile of h{r]) is also 
sketched. Empty circles indicate transitions of second order, while the full ones denote 
transitions of first order. In detail, for 3- and 4-layers, the transitions 30 3R and 
4n 4:R occur continuously by continuously changing the angle 9 between the two in 
plane basis vectors, in analogy to 20 2R (cf. figure Hj), while all other transitions are 
discontinuous. Additionally, by the transitions 3R 3A, 4i? 4A and 5R 5A, 
and by the transitions changing the layer number at ry = 1.53 (3A 20), r] = 10.14 
(2A ^ 30), rj = 30.03 (3A ^ 40), rj = 66.24 (4A ^ 5R) and r] = 123.11 (5A ^ 6R) 
the distance between outermost layers exhibits a certain jump Ah (indicated by thick 
arrows in figure [71). In fact, there is here no continuous transition present between two 
unequal layered phases as in the case of hard spheres * . 

Furthermore, for high densities, the concrete lattice evolves to a continuous such 
that effects due to the concreteness get negligible. This means, electrostatically, that 
each layer of a m-layered structure is completely compensated by a certain part of 
background as much as 1/m of the whole. 

In this paper we have dealt with a system consisting of particles (macroions) 
interacting via the unscreened Coulomb potential and of particles of opposite charge 
(counterions) ^ which are homogeneously smeared out over a hard slit of width L, 
compensating the charge of the macroions. To determine the stability diagram of 
crystalline phases, we have performed lattice sum calculations of a set of candidates. 
As possible candidates we have taken into account phases with up to six layers (m = 
1, • • • , 6) whose primitive cell contains up to eight particles (n = 1, • • • , 8). Additionally, 
we considered the buckling phases from |50j, too. We have analyzed a regime up to 

130 in our investigations. For small densities, we could trace the existence of the 
triangular monolayer lA. Crossing a certain critical density r]c the system buckles and 
evolves to a trilayered structure. This transition density is also calculated analytically 
by applying a Taylor expansion to the lattice sum for small separations. Furthermore 
the evolving of the layer separation from monolayer to trilayer could be characterized as 
h{r]) ^ {t] — rjcY^'^^ qualitatively. Tuning the density upwards, we have noticed diflFerent 
stable bilayered structures, same as Wigner crystals. Beyond the bilayers, we could also 
find stable tri-, four-, five- and six-layers in square, rhombic and triangular bases. The 
final stability sequence for m > 4 reads therefore: mR mA {m + 1)R with a 
remarkable vanish of square-based phases, where the sequence for m = 3 and m = 4 is 

* In the case of bilayered hard spheres, one can achieve a continuous layer increase from 2A to four- 
layered hcp-like and hcp(lOO) phase pH [32l [34j . 
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Figure 7. Order parameter h of all stable crystalline phases. Empty circles denote a 
continuous transition, while the full circles mark a discontinuous one. The transitions 
between different layer numbers, rendered as dashed lines, are also recorded as a first 
order transition except lA 3A. Apart of that, the underlined /i-numbers give the 
limit /i- value (77 00), for the case that no more phase transition to a higher layered 
structure occurs. The dotted line indicates a scale change in 77-axis. 



mD mR mA ^ (m + !)□. While the stability domain of evenly layered phases 
gets larger with increasing m, the stability domain of square phases (□) decreases for 
m > 2 and disappears finally for m > 5. On the other hand the stability domain of 
rhombic (R) and triangular (A) phases increases both with growing m > 2. 

Apart of that, the transitions involved here are all of second order except mR mA 
and mA (m + l)n. The latter takes place discontinuously due to the order parameter 
9 and particle positions (as in the case of nR nA) as well as with respect to h (cf. [7]). 

5. Conclusions 

To summarize: For slit-confined ions in a smeared background, we have determined 
the ground state crystalline lattice as a function of the ion density up to the six-layer 
regime. A complex cascade with buckled, squared and triangular bi-, tri-, tetra-, penta- 
and hexalayers was found. The results are verifiable in systems with classical ions in a 
background including charged colloids, dusty plasmas and classical ions in a trap. One 
important conclusion is that the details of multilayered structures depend crucially on 
the particle-background interaction. More future work is needed to include wall charges, 
wall particle attractions and effects of finite temperature [54] . A detailed understanding 
of the stable crystalline structure as originating from the wall properties is desirable to 
construct filter devices [45] or optical band-gap crystals [55]. 
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Appendix 

The total interaction energy per unit cell of a crystalline unscreened Coulomb system 
can be written as 



where the unit cell consists of n particles of charge q located at r^. The self energy in 
equation f lA.ip stems from the interaction between a particle of the unit cell and its own 
periodically repeated images. The term [/^ in equation f lA.ip . is due to the interaction 
between a particle of the unit cell and all other remaining n — 1 particles of the cell 
including their own images. The convergence involved in these sums is guaranteed by 
the inclusion of a surface neutralizing background for each layer. Following the route of 
Brodka and Grzybowsky (see equations (16a), (166) and (17) of reference [48]), [/^ and 
Uq are given below. Therefore reads 



with 7e = 0.577215665 denoting the Euler-Mascheroni constant, Ko{x) the modified 
Bessel function of the second kind [56] and a^^, bx and by the corresponding x- and 
^/-components of the lattice vectors a and b. Using the components Xij = Xi — Xj, 
Vij = Vi — Vj and Zij = Zi — Zj of the relative separation vector r^j between cell particles 
i and j, can be written as 



(A.l) 




(A.2) 




j>i 
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(A.4) 



for (^/^j, Zij) = (0, 0), where is the digamma function [56]. 

Being interested in the transition from mono- to trilayers, we take as input 
the structure characteristics of the triangular phase lA into the lattice sums f lA.2p - 
(|A3): = 7r/3, b^/a^ = 0.5, by/a^ = \/3/2, 7 = 1, Xu/a^ = 0.5 = Xsa/a^, 
yu/by = 1/3 = ^is/^^x = 1, yis/by = 2/3, p = N/A = ^ = |4 and therefore 

= = ^^^^ • Here we consider for lA a multicell (n = 3) consisting of three 
primitive cells, containing each 1 particle. Thus, for a given 77, the energy function Uc 
depends now only on Zu = hL = ^23- Taking this feature into account, the self energy 
and the cross energy finally read 

U'c = I 4 ^ cos (vrfcrn) Ko (nkmVs^ + 7e - In I (A.5) 



and 



1 3 3 2 

1 y^y^g 
In I Z^Z^ e 



i=l j=l 

j>l 



00 



X 



cos I 27rfc^a±^ 1 K„ I 27rt 



,,A:=1 

.cos(27rA;^^^^^lA'n(27rA; 



1/2 



A 



nl/2 



+4E' 



fc=l 

In 



Xj, 



COS 27r/i;^ /^n 2T^k 



2 n 1/2^ 



a, 



cosh (27r|(/)^j|/i) — cos ( 27r— ^ 



- In2 



(A.6) 
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where = V(^i2 ± &y"i)V«x = ^371(1/3 ±m) = A|, Xf^ = y/ {yis ± hymf/al 



V^(2/3 ± m), A2 = L/a, = P23, P13 = ^L/a,, = J2r]/3V^ = (P23 and 



013 = 2r]/3\^. Before expanding the energy function at /j = 0, we first define 

f{h)^ = Ko (^27rk [a±' + /?2/i2] '^'^ , (A.7) 

where the first four derivatives of /(/i) at /i = are given as follows: 

/(0)± =Ko{27rkX^), (A.8) 
/'(O)^ =0, (A.9) 

nor =-K,{2nkX^)^, (A.IO) 

r (0)± = 0, (A.ll) 

r\0)^ = [Ko {27rkX^) 27rkX^ + 2K, {27rkX^)] (A.12) 

Here, Ki{x) is a modified Bessel function of the second kind [56], too. Using a Taylor 
series and (IA.8[) - (IA.12[) . we now expand Uc{h) from (lA.lj) at /i = and achieve the final 
form of the energy: 

euih) 4 

= -1.960516 - 3.59066877/1^ + 4.968827ry2/i^ +o^^^V^- (A.13) 
rv^ ^ V ' 3 

The last term stems from ([5j), due to interactions with the background, respectively. 
The coefficient —1.960516 corresponds to the static energy per particle of the triangular 
lattice, which was already calculated in [57]. 
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